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ABSTRACT 


An  analysis  is  presented  of  MHD  channel  flow  with  Hall  effect 
and  variable  fluid  properties.  Because  of  the  highly  nonlinear  nature 
of  the  governing  conservation  equations,  a  straight  numerical  analysis 
of  the  equations  is  followed  based  on  the  implicit  finite  difference 
method.  The  analysis  considers  the  two-dimensional  flow  in  the  elec¬ 
tric  field  plane  in  which  the  channel  walls  are  electrode  surfaces. 
Because  of  the  presence  of  the  Hall  effect,  a  transverse  pressure 
gradient  is  induced  in  the  flow  which  necessitates  the  use  of  two  equa¬ 
tions  of  motion  in  the  longitudinal  and  transverse  directions.  The 
analysis  model  employed  is  the  so-called  slender  channel  model  in 
which  the  boundary  layer  forms  of  the  equations  of  motion  are  applied 
across  the  entire  channel  width.  Because  the  boundary  layer  equa¬ 
tions  are  parabolic,  it  is  possible  to  integrate  the  equations  in  a 
forward  marching  manner  by  starting  with  initially  prescribed  veloc¬ 
ity  and  enthalpy  profiles  at  the  channel  entrance.  Numerical  results 
are  presented  for  the  velocity  and  enthalpy  profile  development  from 
initially  assumed  profiles  of  velocity  and  enthalpy  at  the  channel  en¬ 
trance.  Initial  velocity  and  enthalpy  profiles  which  are  typical  of 
both  laminar  and  turbulent  flows  are  considered.  In  addition  to  the 
calculation  of  profiles,  the  voltage  drop  across  electrodes  is  cal¬ 
culated  by  integration  of  the  conductivity  profile  across  the  channel. 
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SECTION  I 

INTRODUCTION 


Two  promising  areas  of  engineering  applications  of 
magnetohydrodynamics  (MHD)  have  been  defined  in  recent  years, 
namely  MHD  power  generation  and  MHD  acceleration.  The  accel¬ 
erator  application  has  received  serious  study  for  several 
years  as  a  technique  for  augmenting  the  performance  of  aero¬ 
dynamic  test  facilities.  The  principal  advantages  of  an 
MHD  augmented  test  facility  such  as  a  shock  tunnel  are: 

(1)  energy  is  added  to  a  moving  stream  resulting  in  test 
section  flows  with  stagnation  pressure  and  temperature 
which  would  exceed  those  achievable  in  a  reservoir  because 
of  material  limitations,  and  (2)  the  energy  is  added  as 
directed  kinetic  energy  rather  than  thermal  energy. 

In  order  to  design  MHD  augmented  accelerators  with 
specified  performance  characteristics,  the  prediction  of  the 
gas  conditions  in  the  accelerating  and  test  regions  is  essen¬ 
tial.  Because  of  the  extremely  complex  physical  phenomena 
occurring  in  MHD  devices,  it  is  safe  to  say  that  at  this 
stage  of  development,  a  complete  theoretical  description  of 
the  flow  is  not  possible.  Because  of  questions  raised  re¬ 
cently  concerning  the  role  of  instabilities  on  the  perfor¬ 
mance  of  MHD  devices,  it  is  at  best  an  open  question  as  to 
whether  the  traditional  analytical  tools  of  the  fluid  mech¬ 
anician,  namely  the  continuum  partial  differential  conser¬ 
vation  equations  coupled  with  the  electromagnetic  field 
equations,  are  adequate  to  describe  the  flow  in  principle. 

In  addition  to  the  question  of  instabilities,  other  physical 
phenomena  which  make  the  analysis  of  MHD  flows  difficult  are 
(1)  the  variation  of  thermodynamic  and  transport  properties 
resulting  in  a  strong  coupling  between  the  fluid  mechanics 
and  electromagnetic  equations,  (2)  the  two  dimensionality  of 
most  flows  of  practical  interest  (small  aspect  ratio  internal 
channel  flows  are  actually  three-dimensional),  (3)  the  high 
nonlinearity  of  the  governing  equations,  and  (4)  the 
possibility  of  turbulence  occurring  for  sufficiently  high 
Reynolds  numbers. 

Any  analytical  investigation  of  such  complex  physical 
phenomena  as  MHD  flows  requires  the  adoption  of  models  which 
do  not  include  all  of  the  physical  phenomena  involved.  For 
an  analytical  model  to  be  useful,  it  must  retain  the  features 
of  those  phenomena  one  wishes  to  focus  attention  on.  For 
example,  a  knowledge  of  the  rate  of  growth  of  wall  boundary 
layers  in  MHD  acceleration  devices  is  important  in  evaluating 
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the  performance  of  the  device  as  well  as  determining  the  ex¬ 
tent  of  the  flow  with  reasonably  uniform  properties  which 
could  be  used  as  a  wind  tunnel  test  flow. 

The  boundary  layer  equations  themselves  can  be  considered 
as  an  analytical  model  which  focuses  attention  on  that  region 
of  the  flow  in  which  the  inertia  and  viscous  terms  are  of 
equal  importance.  Given  the  boundary  layer  equations,  there 
are  then  available  several  approaches  for  extracting  the  in¬ 
formation  on  the  boundary  layer  flow  region  contained  in  the 
equations.  These  methods  include  similarity  solutions,  in¬ 
tegral  methods,  finite  difference  methods,  and  linearization 
methods  of  the  Oseen  approximation  type.  Still  another  method 
which  has  recently  been  applied  by  the  authors  (Ref.  1)  to 
MHD  boundary  layer  flows  in  internal  channels  is  the  slender 
channel  model.  In  this  method,  the  internal  channel  flow  is 
not  divided  into  a  boundary  layer  -  core  flow  region  in 
which  different  equations  are  to  be  applied.  Instead,  the 
boundary  layer  form  of  the  equations  of  motion  are  applied 
across  the  entire  channel  width  with  boundary  conditions 
being  applied  at  the  channel  walls.  This  method  was  also 
recently  applied  by  Kitowski  (Ref.  2)  to  analyze  MHD  flow 
in  the  entrance  region  of  a  parallel  plate  channel. 

In  the  present  investigation,  the  slender  channel  model 
is  used  to  analyze  MHD  flow  in  the  E-plane  (the  plane  con¬ 
taining  the  electric  field  vector)  of  a  two-dimensional 
channel  including  Hall  effect.  The  channel  walls  thus 
correspond  to  the  electrode  walls  of  an  MHD  generator  or 
accelerator.  The  emphasis  Is  placed  on  the  acceleration 
mode  of  operation  in  the  present  study. 
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SECTION  II 

PREVIOUS  WORK 


The  interaction  of  an  electrically  conducting  fluid 
media  with  applied  electric  and  magnetic  fields  has  been 
the  subject  of  intense  study  over  the  past  several  years. 

A  brief  summary  will  be  given  of  those  investigations 
most  pertinent  to  the  present  study. 

The  original  theoretical  investigation  of  magneto¬ 
hydrodynamic  channel  flow  was  conducted  by  Hartmann  (3). 

This  study  was  concerned  with  the  fully  developed  flow 
of  an  incompressible  fluid  between  parallel  plane  walls 
with  a  uniform  magnetic  field  applied  perpendicular  to 
the  flow.  This  is,  of  course,  the  basic  Faraday  generator 
or  pump  geometry  and  a  number  of  authors  have  used  this 
model  to  study  in  detail  the  effect  of  various  modifications 
of  the  Hartmann  problem.  Eraslan  (4)  and  Young  (5)  have 
examined  the  effect  of  variable  electrical  profiles. 

Snyder  (6)  used  an  Oseen  linearization  technique  to  ex¬ 
plore  the  development  of  the  flow  in  the  entrance  region 
of  a  parallel  plate  channel. 

The  equations  for  compressible,  inviscid  MHD  channel 
flow  were  formulated  by  Resler  and  Sears  (Ref.  7)  using  a 
quasi-one-dimensional  approximation.  Since  the  quasi-one- 
dimensional  model  Involves  more  unknowns  than  equations, 
the  streamwise  variation  of  some  of  the  independent  varia¬ 
bles  must  be  specified  before  a  solution  can  be  obtained. 
Special  cases  for  prescribed  variations  of  electric  field, 
magnetic  field,  channel  area,  etc.  have  been  calculated. 

The  results  of  the  quasi-one-dimensional  analysis 
of  compressible  MHD  channel  flow  provides  free  stream 
conditions  that  can  be  used  in  an  analysis  of  the  viscous 
layer  adjacent  to  the  channel  walls. 

The  boundary  layer  flow  adjacent  to  the  electrodes 
in  an  MHD  accelerator  was  studied  by  Kerrebrock  (8). 

This  investigation  was  concerned  with  compressible  boundary 
layer  flow  and  the  effect  of  variable  fluid  properties  was 
included.  Similar  solutions  were  obtained  for  a  certain 
class  of  accelerator  flows  with  constant  free  stream 
temperature.  The  solution  showed  that  the  heat  transfer 
rates  were  much  greater  than  would  be  the  case  for  non-MHD 
flow. 
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The  Slender  channel  model  differs  from  the  usual 
boundary  layer  approach  In  that  the  boundary  layer  equations 
are  applied  across  the  entire  channel  rather  than  Just  in 
the  vicinity  of  the  wall.  In  this  method  the  streamwise 
variation  of  the  velocity,  pressure,  and  temperature  are 
given  by  the  solution  of  the  equations  rather  than  by  a 
separate  inviscld  analysis.  This  technique  was  first 
applied  by  Williams  (9)  to  ordinary  compressible,  viscous 
flow  in  two-dimensional  and  axisymmetric  slender  channels. 

The  method  was  extended  by  Sonnerup  (10)  to  MHD 
flows.  Sonnerup  investigated  the  two-dimensional,  fully 
developed  flow  of  a  compressible,  electrically  conducting 
fluid  in  a  slowly  diverging  channel.  The  confining 
boundaries  of  the  channel  were  the  electrode  surfaces  of 
an  MHD  generator.  By  a  proper  choice  of  channel  geometry, 
the  electromagnetic  forces  and  the  pressure  forces  were 
made  to  cancel.  This,  in  turn,  permitted  the  momentum 
equation  to  be  uncoupled  from  the  energy  equation  and  re¬ 
duced  to  an  ordinary  differential  equation.  An  exact 
solution  to  this  equation  was  obtained  for  the  self- 
slmilar  velocity  profiles.  The  corresponding  temperature 
profiles  were  generally  non-similar. 

Kubin  and  Smolin  (11)  investigated  the  two-dimensional, 
Isothermal  flow  of  an  electrically  conducting  fluid  in 
the  presence  of  crossed  electric  and  magnetic  fields. 

The  form  of  the  governing  equations  was  very  similar  to 
those  of  Williams  (9) .  By  tailoring  the  electric  and 
magnetic  field  distributions,  self  similar  velocity  solu¬ 
tions  were  obtained. 

The  velocity  and  temperature  distribution  for  magneto- 
gasdynamic  acceleration  flow  was  examined  by  Snyder  and 
Maus  (1) .  The  governing  equations  were  formulated  in  polar 
coordinates  for  steady,  two-dimensional,  laminar  flow 
between  slightly  diverging  insulator  walls.  The  assumption 
was  made  that  the  dependent  variables  could  be  expressed 
in  a  separable  form.  This  allowed  the  governing  partial 
differential  equations  to  be  reduced  to  a  pair  of  coupled, 
non-linear,  ordinary  differential  equations.  These 
equations  were  then  integrated  numerically  to  yield  the 
velocity  and  temperature  profiles.  For  certain  ranges  of 
the  electromagnetic  and  gas  dynamic  parameters,  the 
temperature  distribution  exhibited  a  marked  excess  near 
the  wall.  This  was  caused  by  the  large  Joule  heating  and 
viscous  dissipation  effects  in  the  region. 

Kitowski  (2)  Investigated  the  two-dimensional,  com¬ 
pressible  MHD  flow  in  the  entrance  region  of  a  parallel 
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plate  channel.  The  slender  channel  model  was  used  and  the 
boundary  layer  equations  were  written  in  an  implicit  finite 
difference  form  developed  by  Flugge-Lotz  and  Blottner  (12). 
The  analysis  was  concerned  with  the  flow  between  insulating 
walls  and  Included  the  variation  of  transport  properties. 

The  finite  difference  equations  were  solved  numerically  to 
yield  velocity  and  static  enthalpy  profiles  along  the 
channel. 

None  of  the  studies  discussed  to  this  point  have  in¬ 
cluded  the  Hall  effect  where  an  electrical  current  component 
is  induced  in  the  streamwise  direction  due  to  the  electron 
velocity.  This  Hall  current  results  in  a  much  tighter 
coupling  of  the  basic  governing  equations  and  thus  com¬ 
plicates  considerably  the  math^atlcal  problem  associated 
with  the  MHD  flow. 

The  previous  investigations  on  MHD  flow  Including  Hall 
effects  can  generally  be  classified  into  two  categories 
according  to  whether  the  main  emphasis  of  the  study  is  on 
the  electrodynamics  or  the  fluid  mechanics.  In  the  former 
approach  the  gas  dynamical  variables  are  considered 
apecified  and  solutions  are  sought  for  the  electric  po¬ 
tential  and  current  density  distributions.  These  variables 
are  governed  by  elliptic  partial  differential  equations 
derivable  from  Maxwell's  equations  and  the  generalized 
Ohm's  law.  On  the  other  hand,  in  the  latter  approach, 
solutions  are  sought  for  the  fluid  mechanical  equations 
with  the  electrodynamic  variables  specified.  It  is  the 
latter  approach  that  is  taken  in  the  present  investigation. 

The  studies  of  Celinski  and  Fischer  (13)  and  Oliver 
and  Mitchner  (14)  are  examples  of  analyses  where  the 
electrodynamics  of  the  MHD  flow  is  emphasized.  Celinski 
and  Fischer  examined  the  effect  of  finite  size  electrodes 
on  the  distribution  of  electrical  parameters  in  MHD  gene¬ 
rators  for  a  wide  range  of  Hall  parameter  odt.  The  gas 
dynamical  quantities,  velocity  and  temperature  were  assumed 
constant  throughout  the  channel.  Numerical  solutions  were 
obtained  to  the  partial  differential  equations  describing 
the  electric  potential  and  current  density  distributions. 

The  analysis  of  Oliver  and  Mitchner  is  somewhat  similar 
to  that  of  Celinski  and  Fischer  except  that  the  variation 
of  the  electrical  conductivity  was  Included. 

Thw  work  of  Sherman  and  Sutton  (15)  was  one  of  the 
earliest  fluid  dynamic  studies  of  MHD  channel  flow  in¬ 
cluding  Hall  effects.  This  investigation  dealt  with  the 
fully  developed,  incompressible  flow  between  parallel 
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insulating  wall,  i.e.,  Hartmann  flow  with  Hall  effects. 

The  results  of  this  study  showed  that  the  streamwise 
velocity  distribution  changed  from  the  characteristic 
Hartmann  profile  to  a  Ppissuelle-like  profile  with  in¬ 
creasing  Hall  parameter. 

Snyder  and  Eraslan  (16)  studied  the  MHD  flow  of  an 
incompressible  fluid  in  the  entrance  region  of  a  parallel 
plate  channel  including  Hall  effects.  This  analysis  was 
an  extension  of  Ref.  6  and  used  the  same  linearizing 
technique  employed  in  that  study.  The  results  of  this  in¬ 
vestigation  showed  that  the  entrance  length  Increased  with 
increasing  values  of  the  Hall  parameter,  cor. 

Hale  and  Kerrebrock  (17)  were  the  first  to  consider 
both  compressibility  and  the  Hall  effect  in  treating  the 
boundary  layer  on  the  insulating  wall  of  an  MHD  accelerator. 
Two  models  for  the  electrical  conductivity  variation  were 
employed:  an  equilibrium  model  where  the  degree  of  ionir 

zation  is  determined  by  the  local  gas  temperature  and 
pressure,  and  a  non-equilibrium  model  in  which  the  elec¬ 
trical  conductivity  is  strongly  coupled  to  the  electric 
field  strength. 

It  was  found  that  similar  solutions  could  be  obtained 
only  for  very  restrictive  conditions  and  local  similarity 
was  assumed.  The  partial  differential  equations  of  motion 
were  thus  reduced  to  ordinary  differential  equations. 

Some  numerical  difficulty  was  encountered  in  dealing  with 
the  resulting  two  point  boundary  value  problem  and  for 
certain  cases  satisfactory  convergence  at  the  edge  of  the 
boundary  layer  was  not  obtained. 

The  results  obtained  by  Hale  and  Kerrebrock  indicated 
that  the  general  effect  of  Increasing  the  Hall  parameter 
was  to  increase  the  average  current  near  the  wall  leading 
to  a  thinner  boundary  layer,  increased  heat  transfer, 
shear  stress,  and  electrical  losses. 

Hunter  (18)  in  a  recent  study,  analyzed  the  boundary 
layer  on  the  insulating  walls  of  an  MHD  generator  including 
Hall  effects.  The  governing  equations  were  written  in 
modified  Levy-Lees  variables  and  non-similar  solutions 
were  obtained  by  numerically  integrating  the  equations 
using  an  implicit  finite  difference  method.  The  effect 
of  the  boundary  layer  growth  on  the  inviscid  flow  field 
was  accounted  for  through  a  displacement  thickness  correction. 

The  results  of  this  investigation  showed  that  the  Hall 
effect  has  a  rather  pronounced  influence  on  the  flow  in  MHD 
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generators.  It  was  found  that  the  Hall  current  tends  to 
cause  boundary  layer  separation  further  upstream  for  the 
decelerating  flow  in  a  constant  area  generator.  The  energy 
dissipation  due  to  joule  heating  was  underpredicted  when 
the  Hall  effect  was  not  included. 
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SECTION  III 

FORMULATION  OF  THE  PROBLEM 
3 . 1  GOVERNING  EQUATIONS 


The  geometry  for  the  two-dimensional  MHD  flow  in  the 
E-plane  is  shown  in  Figure  1.  The  channel  dimension  in  the 
z-direction  is  assumed  to  be  infinite  thus  reducing  the 
problem  to  a  two-dimensional  flow.  The  channel  height,  W(x), 
is  assumed  to  be  a  specified  function  of  x,  and  the  magnetic 
field,  B^,  is  constant  and  is  in  the  z-direction. 

The  assumption  of  small  magnetic  Reynolds  number  is  made 
which  implies  that  the  equations  of  motion  and  the  magnetic 
field  equation  are  uncoupled.  Thus  the  induced  magnetic 
field  effect  on  the  equations  of  motion  is  neglected.  Using 
the  slender  channel  assumption  means  that  the  complete  equa¬ 
tions  of  motion,  which  are  elliptic,  are  reduced  to  parabolic, 
boundary  layer  form.  This  means  that  the  equations  may  be 
integrated  in  a  marching  fashion  by  starting  with  initial 
profiles. 

With  the  inclusion  of  the  Hall  effect  in  the  equations 
of  motion,  the  transverse  equation  of  motion  must  be  in¬ 
cluded  because  of  the  Hall  current  induced  transverse 
pressure  gradient.  The  governing  equations  may  be  written 
as  follows. 


^p* 

+ 

(n* 

+  j?bJ 

(1) 

Sx* 

by* 

5y* 

ay*  y  ^ 

Sv*  * 

^p* 

(n* 

— )  -  j!b* 

Sx* 

by* 

Sy* 

by* 

ay* 

(2) 

V* 

=  u* 

bP*  ^ 

ax» 

ay* 

bx* 

by* 

Pr  by* 

^  T  *2  ,  T  *2 

3y»  o* 
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^(o*u*)  ^  ^(o*v*)  ^  0 
3x*  ^y* 


Equations  (l)-(4)  are  respectively  the  two  equations  of 
motion,  the  energy  equation,  the  continuity  equation.  These 
conservation  equations  must  be  supplemented  with  an  equation 
of  state.  Ohm's  law  for  the  electric  current,  and  a  speci¬ 
fication  of  viscosity,  p,*,  conductivity,  a*.  Hall  parameter, 
cjdt,  and  Prandtl  number,  Pr,  as  a  function  of  pressure  and 
temperature.  In  the  present  analysis,  the  Prandtl  number 
and  the  heat  capacity  ratio,  7,  will  be  taken  as  constant. 

An  ideal  gas  equation  of  state  will  be  assumed. 

The  equation  of  state  and  the  components  of  electric 
current  from  Ohm's  law  can  be  written  as 


p*  =  (Ili)p*h* 
7 


(5) 


(6) 


(7) 


The  electric  current  components  given  by  Eqns.  (6)  and  (7) 
introduce  the  electric  field  components  as  additional  un¬ 
knowns.  An  equation  for  the  electric  field  potential  can 
be  derived  from  the  relations 


E  *  = 

^  dx* 


E  *  =  ^ 

y  Sy* 


bJ*  dJ  * 

A 


dx 


(8) 


(9) 
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Equation  (9)  is  the  conservation  of  electric  charge  equation 
(current  continuity).  Combining  Eqns.  (6)-(9)  gives 

^  ^  ^  ^^o*  ^  el 

^x*2  ^y*2  3x*  ay*  cr*  Sx*  a*  ay* 

^  ^ 

ay*  ax*  aj  ay*  aj  ^x* 

+  B  *[-^  (v*  +  ojTu*)  -i — ^  (onrv*  -  u*)] 

ax*  ay* 


+  B  *(v*  +  OJTU*) 

z 

+  B^*(£DTV*  -  u*) 


(10) 


where 


1  +  (cur) 


Equation  (10)  is  an  elliptic  equation  describing  the  electric 
field  potential  and  use  of  this  equation  requires  the  speci¬ 
fication  of  boundary  conditions  for  <I)*  or  its  normal  deri¬ 
vative  on  a  closed  boundary.  The  terms  on  the  right  side  of 
Eqn.  (10)  involve  the  fluid  properties  and  thus  there  is  a 
strong  coupling  of  the  electric  field  equations  and  the  fluid 
mechanics  equations.  It  is  seen  that  the  right  side  of  Eqn. 
(10)  becomes  zero  only  if  all  fluid  properties  including 
velocity  are  constant.  In  this  case,  the  equation  reduces 
to  Laplace's  equation  which  is  still  elliptic.  In  order  to 
solve  the  fluid  mechanics  and  electric  field  equations  simul¬ 
taneously,  an  Iteration  procedure  must  be  used  in  which  a 
tentative  solution  to  the  fluid  mechanics  equations  is  ob¬ 
tained  based  on  an  assumed  electric  field  distribution. 

From  this  solution,  the  fluid  property  dependent  terms  in 
Eqn.  (10)  may  be  evaluated  and  a  solution  for  the  electric 
field  potential  obtained.  With  this  corrected  solution  to 
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the  electric  field,  the  fluid  mechanics  equations  can  then 
be  solved  again  with  the  process  being  repeated  as  many 
times  as  desired. 

Such  an  iteration  procedure  would  be  very  time  con¬ 
suming,  and  was  not  used  in  the  present  investigation.  An 
approximate  method  was  used  which  will  be  discussed  in  detail 
later. 

To  account  for  the  divergence  of  the  channel  walls  as 
Illustrated  in  Figure  1,  a  dimensionless  transverse  coordi¬ 
nate  T]  is  defined  as 


T]  = 


y* 

w*(x*) 


(11) 


Transforming  from  (x+,y*)  coordinates  to  (x+,7])  coordinates 
can  be  done  by  utilizing  the  partial  differential  operators 


_i_)  =  -i-)  -  JL  1.) 

Sx*  y*  ^x*  T)  w*  dx*  X* 

(12) 

Sy*  X*  w’^  dri  x* 


Dimensionless  variables  are  defined  as  follows. 


=  si 


X  = 


X 


w 


u  = 


u" 


u. 


V  = 


w  = 


u 


w 


w 


M- 


=  ii!  h  = 


u 


Re  = 


n  u  w 
Po  o  o 


p  = 


p  u 
^o  o 


(13) 


e: 


X 


e: 


u 

o  z 


u  B* 
o  z 


a  = 


H  =  w  B*  -2 
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The  quantities  subscripted  with  o  are  constant  reference 
values  to  be  specified  later.  Rewriting  the  governing 
equations  in  dimensionless  form  and  introducing  the  trans¬ 
formations  of  Eqn.  (12)  gives 

-  +  -A  1.  (p, 

dx  w  dx  ^T]  w  St]  ^x  Rew^  Bt] 


+  - - - K-  +  OTTft  )  +  (COTV  -  u)] 

Re  1  +  (oyr)^  ^  * 


(14) 


pu[^  -  n^^]+£vav^_l^^  1 

^x  w  dx  w  Sti  w  Rew^ 


pu[ 


3x 


-  as - 2 -  |-(*  _  yyr*  )  +  iv  +  oyru)] 

Re  1  -t-  (cot)^  *  y 

-aaSL^]  +  £Z^=u^-t-  (-^  ^)  + 

w  dx  ^ri  w  ^Tj  ^x  Rew"^  Pr  ^t| 


(15) 


Ha^  a 


Re[l  +  (cm-)^] 


I 

{[(4'y  -  CDT$  )  -f  (v  -p  oyru)]' 


(16) 


3pu  _  n  ^  ^pu  ^  ^  3pv 
Sx  w  dx  w  Sti 


(17) 


P  =  (^)ph 
7 


(18) 


Equations  (14) -(18)  are  a  set  of  five  coupled,  nonlinear 
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equations  for  the  five  unknown  quantities  P,  p,  h,  u,  v. 
The  equations  are  parabolic  which  means  that  numerical  in¬ 
tegration  can  be  carried  out  in  a  forward  marching  manner 
by  starting  with  initially  specified  profiles  at  a  given 
X  position  in  the  channel.  Suitable  boundary  conditions, 
to  be  described  in  the  next  section,  must  be  prescribed. 


3.2  BOUNDARY  CONDITIONS 


The  range  of  the  dimensionless  transverse  coordinate 
T]  is  0  ^  T]  ^  1.  Thus  boundary  conditions  must  be  specified 
along  T]  =  0  and  ti  =  1.  The  velocity  boundary  conditions  are 
the  usual  nonslip  conditions.  The  boundary  conditions  on  h 
can  be  specified  either  as  constant  heat  flux  or  constant 
enthalpy.  The  latter  conditions  were  used  in  the  present 
study.  Thus  the  velocity  and  enthalpy  boundary  conditions 
may  be  written  as  follows. 


T^  =  0  u  =  0  v  =  0  h  =  h^j  =  constant 

(19) 

’^  =  1  u  =  0  v  =  0  h  =  h^2  =  constant 

The  electrical  boundary  conditions  are  more  difficult 
to  handle  than  the  velocity  and  enthalpy  conditions.  As  des¬ 
cribed  earlier.  Equation  (10)  which  describes  the  electric 
field  potential  is  elliptic  which  would  require  a  specifica¬ 
tion  of  boundary  conditions  on  a  closed  boundary.  Due  to 
the  strong  coupling  between  the  fluid  mechanics  and  elec¬ 
trical  equations,  a  complex  and' time  consuming  iteration 
procedure  would  have  to  be  used  to  solve  the  coupled  equations. 
An  approximate  but  physically  reasonable  method  of  handling 

the  electric  field  distribution  was  used  in  the  present  in¬ 

vestigation. 

The  total  Hall  current  flowing  in  the  channel  is  given 
by  the  expression 


=  *0*  f 
0 


(20) 


Substituting  for  J*  from  Equation  (6)  and  introducing 
dimensionless  quantities  gives 
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CT  u  w  WB  * 
o  o  o  z 


=  I 


1  +  (out)' 


(1  -  ciyr 


<I> 


X 


1 

+  J  - - - K-  (v  +  ayru)dr| 

Q  1  +  (aiT)^ 


(21) 


At  this  point  two  additional  assumptions  are  introduced. 
The  first  is  that  the  ratio  is  a  constant  which  would 

be  a  valid  condition  for  a  channel  of  diagonal  wall  construc¬ 
tion.  The  second  assumption  is  that  which  is  a  variable 
across  the  channel,  will  be  replaced  by  an  average  value 
which  depends  only  on  x  but  not  t].  Solving  for  from 

Equation  (21)  gives 


4,(x) 


a  u  w  wB^^ 
o  o  o  z 


1 

-  / - - - S-  (v  +  0Cfru)dTi 

Q  1  +  (cut)^ 


J - 2 - -  (1  -  (JUT  ■^)dri 

Q  1  +  (out) 


(22) 


The  total  Hall  current,  itf  and  <I)y/$x  ^^re  taken  as 
constant  parameters  of  the  analysis.  After  profiles  at  each 
station  have  been  calculated,  a  value  of  then  cal¬ 

culated  at  that  station  using  Equation  (22). 


3.3  ANALYSIS  OF  PRESSURE  DISTRIBUTION 


In  applying  the  first  order  boundary  layer  equations 
to  an  external  flow,  the  pressure  distribution  as  given  by 
the  axial  pressure  gradient  dP/dx  is  known  or  at  worst  can 
be  obtained  by  solving  the  inviscid  equations  of  motion  for 
the  given  body  shape  independent  of  the  boundary  layer  equa¬ 
tions.  The  situation  is  more  complex  for  an  internal  channel 
flow  because  the  pressure  distribution  must  be  determined 
concurrently  with  the  solution  to  the  boundary  layer  equations. 
This  is  because  of  the  constraint  of  overall  mass  conservation 
in  which  the  total  mass  flow  rate  is  the  same  at  each  x- 
location.  By  contrast,  the  unconfined  nature  of  an  external 
boundary  layer  flow  means  that  the  inviscid  pressure  distri¬ 
bution  is  not  affected  by  the  growth  of  the  boundary  layer 
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and  consequent  increase  in  boundary  layer  mass  flow. 

A  considerable  effort  was  necessary  to  develop  a  method 
for  calculating  the  transverse  pressure  distribution  at  grid 
column  m+1  in  terms  of  known  variables  at  column  m.  A  first 
approximation  was  made  by  neglecting  the  Inertia  and  viscous 
terms  in  Equation  (15),  the  transverse  equation  of  motion. 
The  body  force  term  was  evaluated  in  terms  of  the  known 
quantities  at  column  m,  and  the  pressure  distribution  was 
calculated  at  m  +  1/2,  the  station  halfway  between  stations 
m  and  m+1.  Thus, 


n 

P(ll)]  =  W  /  f  ]  dT]  +  P(0)] 

m+1/2  ^  m  m+1/2 


(23) 


where  is  the  y  component  of  body  force  evaluated  at 

column  m.  Assuming  the  pressure  at  m  +  1/2  to  be  the 
arithmetic  average  of  the  pressures  at  m  and  m+1  then  gives 


P(n)]m+l  “  ^^^’l^^m+1/2  “  (24) 


The  pressure  at  the  lower  wall,  P(0)]j^^i/2  appearing 
in  Equation  (23)  is  unknown,  however.  An  '  iteration 
must  be  performed  as  follows  to  determine  the  correct  value 
of  P(0)]_,,y«.  An  initial  guess  of  P(0)]„,i yo  is  made  by 
setting  P(0)]m+l/2  =  P(0)]m-  ^ 

This  initial  estimate  of  the  wall  pressure  then  allows  the 
static  pressure  distribution  at  station  m+1  to  be  calcu¬ 
lated  from  Equation  (24).  Knowing  the  pressure,  enthalpy, 
and  velocity  profiles  at  station  m+1  then  allows  the  total 
mass  flow  rate  to  be  calculated.  If  the  calculated  mass 
flow  rate  does  not  agree  with  the  known  mass  flow  rate  based 
on  the  initial  profiles,  an  adjustment  of  the  wall  pressure 
P(0)]m+i/2  is  made  and  the  procedure  repeated  until  con¬ 
vergence  of  the  mass  flow  rate  is  obtained.  In  practice  it 
was  found  that  two  such  iterations  were  sufficient. 

In  carrying  out  the  numerical  analysis,  the  iteration 
was  actually  performed  on  the  unknown  value  of  wall  density 
p(0)]^  ^^2  instead  of  wall  pressure.  This  is  because 
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Equations  (14),  (16),  (18)  were  solved  simultaneously  for 
p,  u,  h  at  station  m+1  and  a  wall  value  of  density  must  be 
assumed  to  initiate  the  numerical  calculation. 


3.4  FINITE  DIFFERENCE  FORMULAE 


Because  of  the  strong  coupling  and  nonlinearity  of 
Equations  (14)-(18),  an  implicit  finite  difference  method 
was  used  to  obtain  an  approximate  solution  to  the  equations 
at  discrete  points  across  the  channel.  For  this  purpose, 
the  two-dimensional  region  in  the  (x,q)  plane  is  divided 
into  a  grid  as  illustrated  in  Figure  2.  The  partial  differ¬ 
entiation  operators  in  the  governing  equations  are  then  re¬ 
placed  by  difference  quotients.  This  allows  the  governing 
equations  to  be  represented  approximately  by  a  set  of  dif¬ 
ference  equations.  Solutions  to  the  difference  equations 
are  obtained  at  intersection  points  of  the  grid  lines.  A 
general  nodal  point  is  represented  by  a  double  subscript 
(m,n) . 


There  are  several  types  of  finite  difference  quotients 
which  may  be  used  to  approximate  the  partial  derivatives. 

Only  implicit  finite  difference  quotients  are  used  in  the 
present  investigation  since  the  stability  of  the  finite 
difference  method  is  enhanced  by  using  implicit  differences. 
With  the  implicit  difference  method,  solutions  along  a  grid 
column  m+1  can  be  obtained  by  simultaneous  solution  of  a 
system  of  linear  algebraic  equations,  utilizing  known  values 
of  the  variable  at  the  previous  m  column.  The  system  of 
equations  is  of  tridiagonal  form  which  means  an  efficient 
algorithm  can  be  utilized  not  involving  the  inversion  of  a 
large  matrix. 

Two  finite  difference  schemes  may  be  used  to  construct 
the  difference  equations:  the  ordinary  implicit  method  and 
the  Crank-Nicolson  method.  The  ordinary  implicit  scheme  is 
developed  by  using  backward  difference  quotients  to  approxi¬ 
mate  partial  derivatives  at  m+l,n.  The  Crank-Nicolson  im¬ 
plicit  scheme  utilizes  central  difference  quotients  to  approx¬ 
imate  partial  derivative  operators  at  m+l/2,n.  The  trun¬ 
cation  error  is  smaller  for  the  Crank-Nicolson  method  than 
for  the  ordinary  implicit  method.  If  f  and  g  represent  two 
arbitrary  dependent  variables,  the  Crank-Nicolson  differences 
may  be  written  as  follows. 


dx 


Vl,n  -  *m,n  1_  ,  ,2 

ax  24 


(25) 
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M 

^T) 


1 

4Ati 


^^m,n+l  ”  ^m+l,n+l  ^m+l,n-l^ 


-  i  +  ••• 


2A'n 


_  _  2f  +  f  -4-  f 

2  ^  m,n+l  m,n  m,n-l  m+l,n+l 


^^m+l,n  '*'  ^m+l,n-l^  ”  g  ^xxrjri 


(26) 


12 


(ATi)^f 


nnnTi 


+  •  •  • 


(M)2  =  ^^m.n+l  ~  ^m.n-1^  ^^m+l,n+l  ~  ^m+l.n-l^ 
2Aii  2A'n 


(27) 


(Ax)‘ 


(f  -  f  f  )  - 

'  XT]  T]  XXT]' 


(Ax)' 


f  f  + 

T]  T]T1T] 


(28) 


M.  ^  =  i  |^^^ro,n+l  ~  ^^m+l,n+l  ~  ^m+l,n-l^ 

St]  ^ri  2  2At]  2A'n 


fe  -  e  )  (f  -  f  ) 

^  ^m,n+l  ^m,n-l  m+l,n+l  m+l,n-l‘^j 

2At]  2A'n 


(Ax) 


(fe  -2fe  +ef  ) 

g  '  T]^XXT]  XT]^XT]  ®T]  XXT]'^ 


(29) 


_  (All) 
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Sx  Ax 


(30) 


f  M 
ari 


m,  n 


1 

4At) 


(g 


m,  n+1 


®m,n-l  ^m+l,n+l 


^m+l,n-l^  +  . . . 

(31) 


The  finite  difference  approximations  of  Equations  (25) 
(31),  when  substituted  into  Equations  (14),  (16),  (18),  re¬ 
sult  in  a  set  of  linear  algebraic  equations  for  the  unknown 
values  of  u,  h,  and  p  at  column  m+1.  The  matrix  of  coeffi¬ 
cients  for  this  set  of  equations  is  the  tridiagonal  form, 
and  the  development  of  the  algorithm  for  solving  the  set  of 
equations  is  contained  in  the  Appendix. 

After  solving  for  the  u,  h,  p,  profiles  at  column  m+1, 
the  pressure  distribution  in  the  transverse  direction  was 
calculated  from  Equation  (25).  Additionally,  the  v  profile 
at  column  m+1  is  calculated  from  Equation  (17),  the  con¬ 
tinuity  equation. 
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SECTION  IV 

DISCUSSION  OF  RESULTS 


A  total  of  17  different  computer  runs  were  made  for 
various  combinations  of  the  parameters.  The  ranges  of  physi¬ 
cal  parameters  for  which  numerical  results  were  obtained  are 
shown  in  Table  I  and  Table  II.  For  each  combination  of 
parameters,  three  sets  of  curves  are  presented.  These  three 
sets  of  curves  show  respectively  (1)  the  centerline  distri¬ 
butions  of  enthalpy,  velocity,  and  pressure,  (2)  velocity  and 
enthalpy  profiles  at  x/Wq  =  4.8,  and  (3)  velocity  and  enthalpy 
profiles  at  x/w^  =  9.8. 

The  Group  I  runs  compare  the  influence  of  initial  profiles 
of  velocity  and  enthalpy  on  the  flow  development.  Initial 
condition  1  corresponds  to  parabolic  profiles  of  velocity  and 
enthalpy  at  inlet  and  initial  condition  2  corresponds  to 
quartic  profiles  of  velocity  and  enthalpy  at  inlet.  The 
quartic  profiles  are  flatter  profiles  than  the  parabolic 
profiles  and  would  simulate  turbulent  initial  profiles. 

For  each  computer  run,  the  calculations  were  done  from 
the  channel  entrance  to  a  downstream  distance  of  approximately 
10  channel  widths.  This  distance  was  sufficient  to  establish 
the  trends  of  the  calculations  and  represents  a  compromise  of 
the  computer  time  that  was  available  between  computing  for  a 
small  number  of  cases  over  a  large  channel  distance  and  a 
large  number  of  cases  over  a  short  channel  distance.  A  few 
runs  were  made  over  a  distance  of  over  100  channel  widths 
with  no  difficulties  encountered  in  the  program. 

It  was  found  that  the  method  of  handling  the  electrical 
boundary  conditions  described  earlier  became  very  time  con¬ 
suming.  Since  the  amount  of  computer  time  available  was 
limited,  all  calculations  were  performed  assuming  constant 
values  of  the  transverse  and  Hall  currents,  J  and 
respectively.  ^ 

The  results  for  the  Group  I  calculations  will  be  dis¬ 
cussed  first.  A  detailed  study  of  the  Group  I  results  con¬ 
tained  in  Figures  3-32  revealed  that  the  following  conclusions 
may  be  drawn. 

1}  The  Hall  current  had  essentially  no  effect  on 

the  centerline  distribution  of  velocity,  enthalpy, 
and  pressure. 
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TABLE  I 


RANGE  OF  PHYSICAL  PARAMETERS  FOR  GROUP  I  RUNS 


Fixed  Values: 

Ha  =  100 

R.  =  0.12  X  10® 

M  =  2.0 

=  0.73 
r 

7  =  1.4 

Jy  =  0.1 

Case 

Jx 

^1 

H^2 

Initial  Conditions* 

1 

0 

0.04 

0.04 

1 

2 

0.05 

0.04 

0.04 

1 

3 

0.10 

0.04 

0.04 

1 

4 

0.05 

0.08 

0.08 

1 

5 

0.05 

0.04 

0.08 

1 

la 

0 

0.04 

0.04 

2 

2a 

0.05 

0.04 

0.04 

2 

3a 

0.10 

0.04 

0.04 

2 

4a 

0.05 

0.08 

0.08 

2 

5a 

0.05 

0.04 

0.08 

2 

♦Initial  condition  1  corresponds  to  parabolic  profiles  of 
velocity  and  enthalpy  at  inlet. 

Initial  condition  2  corresponds  to  quartic  profiles  of 
velocity  and  enthalpy  at  inlet. 
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TABLE  II 

RANGE  OF  PHYSICAL 

PARAMETERS 

FOR  GROUP  II  RUNS 

Fixed  Values: 

Hw2  =  0-04 

=  0.05 

7  =  1.4 

J  =  0.10 

INITIAL  PROFILE  1 


Case 

Ha 

Re 

Pr 

M 

6 

50 

1.2  X  10® 

0.73 

2.0 

7 

100 

1.2  X  105 

0.73 

2.0 

8 

100 

0.6  X  105 

0.73 

2.0 

9 

100 

0.55  X  105 

0.73 

2.0 

10 

100 

1.2  X  105 

0.50 

2.0 

11 

100 

1.2  X  105 

1.00 

2.0 

12 

100 

1.2  X  105 

0.73 

3.0 
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2)  The  shapes  of  the  initial  velocity  profiles  were 
essentially  preserved  over  the  distance  of  10 
channel  widths  .for  which  the  calculations  were 
performed.  Even  though  the  velocity  increased  and 
the  enthalpy  decreased  on  the  centerline,  the 
initially  parabolic  profiles  still  had  a  parabolic 
shape  after  10  channel  widths  and  similarly  for 
the  initially  quart Ic  profiles. 

3)  A  slightly  larger  acceleration  occurs  for  initial 
profiles  which  are  parabolic  compared  to  Initially 
quart ic  profiles.  This  difference  was  slight, 
however,  amounting  to  about  2  to  3  percent  over  a 
distance  of  10  channel  widths. 

4)  The  pressure  drop  increases  slightly  as  the  wall 
enthalpy  increases  for  both  initially  parabolic  and 
quartic  profiles.  The  pressure  drop  is  larger  for 
an  initially  parabolic  profile  than  for  an  initially 
quartic  profile.  Also,  the  centerline  enthalpy  de¬ 
creases  faster  for  an  initially  parabolic  profile 
than  for  a  quartic  initial  profile. 

The  results  for  the  Group  II  calculations  are  contained 
in  Figures  33-53.  A  detailed  study  of  these  curves  indicates 
that  the  following  conclusions  are  valid. 

1)  The  acceleration  on  the  centerline  increases  with 
increasing  Hartmann  number  as  one  would  expect. 

For  example,  with  all  other  conditions  the  same, 
a  Ha  value  of  100  gives  a  14  percent  increase  in 
centerline  velocity  over  10  channel  widths  compared 
to  a  4.5  percent  Increase  for  a  Ha  value  of  50. 

2)  A  smaller  value  of  Reynolds  number  results  in  a 
greater  acceleration  on  the  centerline. 

3)  The  variation  of  Prandtl  number  over  the  range 
0.73  ^  Pr  ^  1.0  had  essentially  no  effect  on  the 
acceleration  characteristics  of  the  flow.  No 
readable  change  in  centerline  values  of  velocity, 
enthalpy,  pressure  can  be  observed  over  this  range 
of  Prandtl  number. 

4)  Increasing  the  Mach  number  results  in  a  decrease 
of  acceleration  on  the  centerline. 
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Applied  B  Field  is  in  the  z-Direction. 


Figure  1.  E-Plane  MHD  Channel  Flow  Geometry 
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Finite  Difference  Grid  Geometry. 


28 


AEDC-TR-70-26 


AEDC-TR-70-26 


AEDC-TR-70-26 


gure  16.  Velocity  and  Enthalpy  profiles  at  x/w^  = 


4.8 
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Figure  20.  Velocity  and  Enthalpy  Profiles  at  x/w^  =  9.8. 
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igure  26.  Velocity  and  Enthalpy 
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Ha  =  100 

Re  =  1.2  X  10® 


Figure  45.  Centerline  Distributions  of  Velocity,  Enthalpy,  Pressure. 
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Figure  48.  Centerline  Distributions  of 
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Figure  50.  Velocity  and  Enthalpy  Profiles  at  x/w^  =  9.8. 
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Figure  51.  Centerline  Distributions  of  Velocity,  Enthalpy,  Pressure. 


AEDC-TR-70-26 


UNCLASSIFIED _ 

DOCUMENT  CONTROL  DATA  -  R  &  D 

(Security  c/<ii///caf/on  of  title,  body  of  mbMtrtel  mnd  Indmxttg  annotation  muat  be  entered  whan  the  o¥erall  repoff  la  cta§mlllmd} 


9.  REPORT  TITLE 

INTERNAL  MHD  CHANNEL  FLOWS  INCLUDING  HALL  EFFECT  AND 
VARIABLE  FLUID  PROPERTIES 

4.  DESCRIPTIVE  NOTES  ^7Vp«  at  report  end  Inclualve  datee) 

Final  Report  November  26,  1969 

8<  AUTHORIS)  (Flret  namet  middle  Inittml,  laet  name) 

William  T.  Snyder  and  James  R.  Mans,  The  University  of  Tennessee 
Space  Institute 


10.  OISTRIBUTION  ITATEMENT 

This  document  has  been  approved  for  public  release  and 
sale;  its  distribution  is  unlimited. 


IS.  ABSTRACT 

An  analysis  is  presented  of  MHD  channel  flow  with  Hall  effect  and 
variable  fluid  properties.  Because  of  the  highly  nonlinear  nature  of 
the  governing  conservation  equations,  a  straight  numerical  analysis  of 
the  equations  is  followed  based  on  the  implicit  finite  difference 
method.  The  analysis  considers  the  two-dimensional  flow  in  the 
electric  field  plane  in  which  the  channel  walls  are  electrode  sur¬ 
faces.  Because  of  the  presence  of  the  Hall  effect,  a  transverse 
pressure  gradient  is  induced  in  the  flow  which  necessitates  the  use 
of  two  equations  of  motion  in  the  longitudinal  and  transverse 
directions.  The  analysis  model  employed  is  the  so-called  slender 
channel  model  in  which  the  boundary  layer  forms  of  the  equations  of 
motion  are  applied  across  the  entire  channel  width.  Because  the 
boundary  layer  equations  are  parabolic,  it  is  possible  to  integrate 
the  equations  in  a  forward  marching  manner  by  staining  with  initially 
prescribed  velocity  and  enthalpy  profiles  at  the  channel  entrance. 


12.  SPONSORING  MILITARY  ACTIVITY 

Arnold  Engineering  Development 
Center,  AFSC,  Arnold  Air  Force 
Station,  Tennessee  37389 


II.  SUPPLEMENTARY  NOTES 

Available  in  DDC. 


6.  REPORT  DATE 

February  1970 


BA.  CONTRACT  OR  GRANT  NO. 

F40600-67-C-0001 ,  Task  I 

b.  PROJECT  NO. 

8950 
Task  12 

d.  Program  Element  6240533F 


7«.  TOTAi.  NO.  OF  PACES  7b.  NO.  OP  REPS 

.88  18 
M.  ORIGINATOR*S  REPORT  NUMBER(S) 

AEDC-TR-70-26 

ftb.  OTHER  REPORT  NO(S)  (Any  othernumbere  that  mmy  be  aeelaned 
thie  teport) 

N/A 


2A.  REPORT  SECURITY  CLASSIFICATION 

UNCLASSIFIED 

zb.  CROUP 


I.  ORICINATINC  ACTIVITY  (Corporalt  Author) 

The  University  of  Tennessee  Space  Institute 
Tullahoma,  Tennessee  37388 


DD  ,''r..l473 


BNCLASSiriED 

Security  Classirication 


UNCLASSIFIED 


curtty  Ciatslfication 


magnetohydrodynamics 
mathematical  analysis 


Anc 

ArwUArtTM 


UNCLASSIFIED 

Security  ClaBaifleatlon 


